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The Burridge-Knopoff model of earthquake faults with viscous friction is equivalent to a van der Pol- 
FitzHugh-Nagumo model for excitable media with elastic coupling. The lubricated creep-slip friction law 
we use in the Burridge-Knopoff model describes the frictional sliding dynamics of a range of real materi- 
als. Low-dimensional structures including synchronized oscillations and propagating fronts are dominant, 
in agreement with the results of laboratory friction experiments. Here we explore the dynamics of fronts in 
elastic excitable media. 



> 

in 

CO 

o 
in 
o 

0\ 1. Introduction 

G\ ■ 



The Burridge-Knopoff model [Burridge & Knopoff, 
,-^1967] mimics the interaction of two plates in a ge- 
i ological fault as a chain of blocks elastically cou- 
^ pled together and to one of the plates, and subject 
^ to a friction force by the surface of the other plate, 
O such that they perform stick-slip motions — Fig. [TJa. 

This simple system reproduces some statistical fea- 
* tures of real earthquakes [ Carlson & Langer, 1989p 
/\ such as the Gutenberg-Richter power-law distribu- 
tion [ Gutenberg & Richter, 1956 ], considered an ex- 
ample self-organized criticality, I Carlson et al., 1994 ] 
which obviously involves a large number of degrees 
of freedom. However, rec ent laboratory experiments 
[ Rubio & Galeano, 1994 ] that attempt to reproduce 
these dynamics in a real stick-slip dynamical system 
consisting of an elastic gel sliding around a metallic 
cylinder, have shown that low-dimensional phenom- 
ena are more robust in reality, and, although proven 
to be unstable in the Burridge-Knopoff model, do 
show up in the laboratory. 

In a different realm, excitable media are usually 
studied using the model of van der Pol, FitzHugh, 
and Nagumo [van der Pol & van der Mark, 1928; 



physicists and mathematicians working in the area 
of nonlinear science because of the apparent univer- 
sality of many features of their complex spatiotem- 
poral properties [ Meron, 1992p. 



We have shown [Cartwright et al., 1997] that a 



Burridge-Knopoff model with a lubricated creep-slip 
friction force law showing viscous properties at both 
the low and high velocity limits (Fig. [lfa) is a type 
of van der Pol-FitzHugh-Nagumo excitable medium 
in which the local interaction is elastic rather than 
diffusive. We have investigated the dynamics of 
the model and have shown that its behaviour is 
dominated by low-dimensional structures, including 
global oscillations and propagating fronts. Here we 
investigate further the dynamics of elastic excitable 
media, and focus on the behaviour of the fronts. 



2. The Model 

Our elastic excitable medium model may be written 
HCartwrightef aZ., 1997p 



dt 2 



c 2 ^-(x 



vt) 



70 l M 



(1) 



FitzHugh, 1960| ; [FitzHugh, 196 1\ Nagumo et al., 
1962 J. This model normally includes only diffusive 
coupling. Originally from physiology and chemistry, 
excitable media have also captured the attention of 



where, in the language of frictional sliding, t) 
represents the time-dependent local longitudinal de- 
formation of the surface of the upper plate in the 
static reference frame of the lower plate, cf>(dx/dt) = 
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Figure 1: (a) The Burridge-Knopoff model, (b) The velocity weakening stick-slip friction law of Carlson & 
Langer [Ff(v) = F a sgn(v)/(l + \v\), where v is the velocity of the block] (solid line) and the Burridge-Knopoff 
type creep-slip friction law we use (dashed line), showing the threshold (vertical line) where the block starts 
to slip. 



(dx/dt) 3 /3 — dx/dt is the friction function, as the 
dashed line in Fig. [l|b, 7 measures the magnitude of 
the friction, c is the longitudinal speed of sound, and 
v represents the pulling velocity or slip rate. Com- 
pare this with the discrete Burridge-Knopoff model 
from whence Eq. ([!]) may be derived in the continuum 
limit 



dt 2 



= kc(xi+i—2xi+Xi-i)—kp(xi—Vt)—Ff 



dxi 
~~dt 



(2) 

where Xi is the departure of block i from its equi- 
librium position. It has been noted [Carlson et al., 
1994] that in some cases discrete Burridge-Knopoff 
models fail to attain a well-defined continuum limit. 
This is not a consequence of any numerical instabil- 
ity in computer simulations, but of the hyperbolic na- 
ture of the equation and of the shape of the friction 
force commonly used. In such cases, Eq. (|j) should be 
considered as a symbolic representation of the well- 
defined discrete dynamics of Eq. (|2|). From Eq. (|l|) 
we can obtain an expression for the local velocity 
V> = dx/dt of the interface that gives us the model 
written as a couple of differential equations of first 
order in time 
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3. Front Dynamics 

The dynamical behaviour of the elastic excitable 
medium mod el has been reported in detail in a pre- 
vious paper QCartwright et al., 1997 ]. It is notable 
that much of the chaotic behaviour shown by the 
discrete versions of the Burridge-Knopoff models 
built on the ba sis of a monotonic velocity weaken- 
ing friction law [ Rice, 1993 ; Schmittbuhl et al., 1993; 



Xu & Knopoff, 1994 ] becomes more organized in our 



case. Global oscillations and propagating fronts dom- 
inate a large proportion of the relevant parameter 



space. Moreover, the global oscillations show inter- 
esting instability mechanisms leading to the appear- 
ance of the propagating fronts. Among these instabil- 
ities, the most interesting is perhaps the occurrence 
of a period-doubling bifurcation at a finite spatial 
wavelength. As a result of this bifurcation, the glob- 
ally synchronized oscillatory medium breaks into a 
finite number of equidistant pacemaker zones from 
which pairs of counterpropagating fronts are emit- 
ted. Fronts emitted from neighbouring pacemakers 
annihilate upon collision and the annihilation point 
then becomes an emitting centre for the next gener- 
ation of fronts; the whole process repeats after two 
iterations. These results are graphically summa- 
rized by Fig. H|, where the transient evolution of a 
nearly uniform initial state is shown for v just above 
the period-doubling instability. In the first stages of 
the evolution the dynamics is dominated by synchro- 
nized global oscillations. The instability then grows, 
giving rise during a certain time interval to a tran- 
sient period-doubled structure. Finally this struc- 
ture decays into a set of propagating fronts (or prop- 
agating pulses, since a pair of neighbouring fronts 
can be considered a pulse). This bifurcation was dis- 
cussed in more detail in a previous work [Cartwright 
et al, 1997]. 

Here we focus on the properties of the propagat- 
ing front regime with special emphasis on the se- 
lection mechanisms for the front velocity and spa- 
tial configuration. We suppose a solution of the type 

ip(x,t) = f(z), where z = x/v + 1, and v is the front 
velocity. This together with the further rescaling 

z — zj y/l — c 2 /v 2 leads to 



-r-j+Ml/ -l)-r + / = f, 

dz z dz 



(5) 



which is the van der Pol equation with the nonlin- 
earity rescaled by fx = 7/ y/1 - c 2 /v 2 . The propagat- 
ing fronts are then periodic solutions of the van der 
Pol equation. The parameter li is undefined until the 
value of the front velocity v is chosen. However, we 
know that the period of the solution is a function T = 
T(fi) of fi: in the limit of large /1, T behaves as T = 
kji + Ojfi- 1 ), where k = 3+ {v 2 - 1) ln[(4- v 2 )/{\ -v 2 )] 
[jFeingold et al., 19880. Since this period should be 




Figure 2: Evolution of global oscillations into propagating fronts via period doubling; a spatiotemporal plot of 
ip(x,t), time is horizontal left to right and space is vertical. 



commensurate with the sys tem size S, we have the 
condition nT(p(v)) — S/ (v-y/l ~ c 2 /v 2 ), where n is an 
integer, to select the allowed front velocities, which 
in the large fi limit gives us the quantizing condition 
v = S/ (nkj). The integer n can be interpreted as the 
total number of pulses propagating in the system. 
Because Eq. (|5|) has bounded solutions only if v 2 > c 2 , 
the propagating fronts are supersonic. Since our an- 
alytical estimations predict that v becomes smaller 
than c and approaches zero as the number of pulses 
in the system is increased, there should be a maxi- 
mum number of pulses allowed in the system. How- 
ever, numerical simulations show that it is possible 
to find solutions composed of an arbitrary number of 
propagating pulses with a proper choice of the ini- 
tial conditions. A question then arises as to what is 
the long term behaviour of these solutions when the 
number of pulses exceeds the maximum allowed by 
the restrictions on the front velocity. There are three 
simple scenarios logically compatible with the anal- 
ysis: 

1. Some of the pulses are annihilated by the dy- 
namics. 

2. Since the configuration cannot propagate 
rigidly at a compatible uniform speed, nonde- 
caying fluctuations of the velocity of individual 
pulses should be observed. 

3. The system evolves into a uniformly propagat- 
ing solution with speed approaching the singu- 
lar limit v = c. In this limit the pulses be- 
come discontinuous, so that Eq. (|]) becomes 
ill-defined. The original system of blocks and 
springs is thus no longer well represented by 
the partial differential equation Eq. (|l|), to 



which our analytic estimations pertain, but the 
discrete effects present in Eq. ^ prevail. 

The first scenario occurs, but only during a tran- 
sient phase: at long times states with arbitrary num- 
bers of pulses are reached and maintained in time. 
Examples of how the second scenario could be real- 
ized will be briefly discussed below. We first focus 
on the third scenario, for which we find abundant 
numerical evidence. Figure ||| summarizes this. In 
upper frame of the first column we see a propagat- 
ing solution composed of only one pulse travelling 
around the system to the right. Using the spatial 
coordinate i as a parameter we have drawn below 
the phase-space diagram of such a solution. This di- 
agram shows geometrically that the spatial depen- 
dence of the solution is accurately described as the 
periodic solution of a van der Pol-FitzHugh-Nagumo 
equation. The second column shows a two-pulse solu- 
tion propagating to the left at a speed approximately 
one half that of the single-pulse solution, and thus 
closer to the sound velocity c, as expected from the 
analytical estimations. The phase-space projection 
below shows general agreement with the prediction 
of Eq. @, although a small discrepancy is already 
visible. Finally, the third column in Fig. || shows a 
breakdown of the continuum model. A train of six 
pulses travelling at constant speed has been gen- 
erated, which is forbidden within the framework of 
the continuum analysis. They travel to the left at 
speed essentially c. Notice the sharp gradients at 
both the leading and trailing edges of the pulses. If 
we refine the numerical spatial discretization, this 
just increases the gradients, so that we may expect 
discontinuities in the continuum limit. The discrete 
model Eq. ^ introduced in the computer should 




Figure 3: Changing the initial conditions one can obtain configurations with different numbers of pulses. In 
the first row we show ip against x, and in the second rj against ip. In the second row the points are data from 
numerical integration, the continuous line is the limit cycle of Eq. (||) with fj, calculated from the parameters 
and the observed v, the dashed line is the nullcline, and the diamond shows the position of the unstable fixed 
point. First column: 1 pulse (2 fronts), velocity v = 4.78 ± 0.64; has the form of the van der Pol— FitzHugh— 
Nagumo limit cycle. Second colum: 2 pulses (4 fronts), velocity v = 2.47 ±0.35; shows small deviations from the 
van der Pol-FitzHugh-Nagumo continuum limit. Third column: 6 pulses (12 fronts), velocity v = 1.10 ± 0.06; 
is clearly outside the continuum limit. The solitary pulse in the first row is travelling to the right; the others 
to the left. Note from the velocity of the solitary pulse that there cannot be five or more pulses in the contin- 
uum limit: since v oc 1/n, they would then have a subsonic velocity. Parameters are 7 = 3,^ = —0.9, c = 1, 
system size L = 60, discretization dx = 60/1024. Integration is by fourth-order Runge-Kutta in time, and the 
second-order finite differences implied by Eq. @ in space. 




Figure 4: The formation of phase fronts between regions of synchronous oscillations; spatiotemporal plots of 
ip(x, t), time is vertical and space is horizontal. 
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better represent the physical system of springs and 
blocks than the partial differential equation Eq. ([!]). 
The phase-space diagram clearly shows a strong 
departure from the solutions of the van der Pol— 
FitzHugh-Nagumo equations. A well-defined phase- 
space structure appears, which should be understood 
as a property of the purely discrete system, different 
in this case from the van der Pol-FitzHugh— Nagumo 
phase-space structure. Discreteness effects have 
been studied i n several models with velocity weak- 
ening friction [Rice, 1993; gchmittbuhl et al., 1993; 



Ku & Knopoff, 1994| ; |Galeano et al., 19980 . 

When v has been taken very far from the slipping 
threshold the behaviour locally is of autonomous re- 
laxation oscillations, and a slightly different kind 
of propagating front dynamics occurs. By setting 
v = we recover the symmetric form of the van 
der Pol equation, supplemented in our case with the 
elastic spatial term. Now a Floquet stability analy- 
sis of the homogeneous oscillatory state shows the 
development of a long wavelength instability as v 
decreases. The evolution of this instability is dis- 
played in Fig. Q One can roughly describe this be- 
haviour as the formation of spatial regions where the 
medium undergoes almost synchronous relaxation 
oscillations. Neighbouring regions, however, oscil- 
late in antiphase and a phase-change front defining 
the border between them travels around the system. 
On the left we can see some of these fronts moving 
at different velocities during a transient stage, while 
on the right is displayed an asymptotic state where 
the fronts have reached an equilibrium configura- 
tion in which the phase distribution travels rigidly. 
This behaviour is reminiscent of the phase dynam- 
i cs found in the complex Ginzburg-Landau eq uation 
HMontagne et al, 1996j ; |Montagne let al, 199?] ]. 

A closer look at Fig. ^| indicates that the spa- 
tial regions that we describe as synchronized, are in 
fact regions where the excitation or relaxation prop- 
agates extremely fast. This is visualized in the pic- 
ture as a very slight tilt of the bands representing 
the oscillation. Notice that the projection of these 
bands shrinks considerably near the places where 
the phase jumps. By taking spacelike slices of these 
pictures we obtain the instantaneous configuration 
of the excitation field. Such configurations look like 
travelling pulses varying wildly in size and velocity. 
This behaviour can be considered an extreme mani- 
festation of scenario 2 above. 



4. Discussion 

A lubricated friction law in the Bu rridge-Knopof f 
model can be j ustified both by the ory [ Persson, 1995 ; 
Persson, 1997|], and experiments iHeslot et al., 1994 ; 
Brechet & Estrin, 1994j; [Kilgore et al., 1993] ; Demirel 
& Granick, 1996; [Budakian et al, 19980. More- 



over, studies of peeling adhesive tape [Hong & Yue, 
1995], of Saffman— Taylor fracture in viscous finger- 



ing [Kurtze & Hong, 1993], and of the Portevin-Le 
Chatelier effect [ [Kubin & Estrin, 1985| ; Lebyodkin 
et al. , 1995] lead to the same form of friction law 
as we use here. Certainly our model well repre- 
sents the qualitative characteristics of the laboratory 
stick-slip dynamics experiments [Rubio & Galeano, 
1994] referred to above, and might have relevance 
to the present debate on shear stress and friction 
in real geolog i cal faults [51eep & Blanpied, 1992 ; 
Melosh, 19961 |Cohen, 1996P - In the context of ex- 
citable systems elastic coupling has been left aside, 
because in the chemical and biological systems stud- 
ied up to now the coupling is diffusive. However, an 
elastic excitable medium can be realized as an active 
transmission line or optical waveguide [Cartwright 
et al., 1997]. 
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